cd "C:\Users\maur3\Dropbox\PC (2)\Desktop\Gate project - commercial agriculture\Productivity paper\dta"
global temp "C:\Users\maur3\Dropbox\PC (2)\Desktop\Gate project - commercial agriculture\Productivity paper\temp"


************************************************************************
*SFA-TRE models (by country) - Tables 5, 6, & A1
************************************************************************
use panel_1013Final.dta, clear

global clust 	"cll"
global idd 		"idd"
global ceth 	"prop10 traderdens ac_ent hirlmarket  cred agriext agritrain area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global cmlw 	"prop10 traderdens ac_ent hirlmarket  cred agriext agritrain area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global cngr 	"prop10 traderdens ac_ent hirlmarket  cred agriext agritrain area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global ctzn 	"prop10 traderdens ac_ent hirlmarket  cred 		   agritrain area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global cugd 	"prop10 traderdens ac_ent hirlmarket  cred         agritrain area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global cghn 	"prop10 traderdens ac_ent hirlmarket  cred agriext 			 area_fertilizer area_nonfarminc lnntl_mean_harm parcels_owned ac_diver NonfIncShare agehead female ac_headeduc hhsize "
global othinp 	"pesticides fertilizer seedmarket hh_any_parcel_irrigated"


****************************************
*ETH
****************************************
{
preserve
keep if country=="ETH"
duplicates tag $idd, gen(tag)
drop if tag==0
tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
sfpanel lnout c.lnA##c.lnA c.lnL##c.lnL c.lnA#c.lnL $othinp tt2 tt3, distribution(exp) vce(robust) usigma($ceth) vsigma() model(tre) marginal itera(100) rescale tolerance(1e-5) nrtolerance(1e-5) 

local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,jlms
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(ETH_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "ETH_tl9.dta", replace
restore
}
*
****************************************
*MLW
****************************************
{
preserve
keep if country=="MLW"
duplicates tag $idd, gen(tag)
drop if tag==0
tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
sfpanel lnout c.lnA##c.lnA c.lnL##c.lnL c.lnA#c.lnL $othinp tt2 tt3 tt4, distribution(exp) vce(robust) usigma($cmlw) vsigma() model(tre) marginal itera(100) rescale tolerance(1e-5) nrtolerance(1e-5) 

local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,jlms
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(MLW_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "MLW_tl9.dta", replace
restore
}
*
****************************************
*NGR
****************************************
{
preserve
keep if country=="NGR"
duplicates tag $idd, gen(tag)
drop if tag==0
tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
sfpanel lnout c.lnA##c.lnA c.lnL##c.lnL c.lnA#c.lnL $othinp tt2 tt3 tt4, distribution(exp) vce(robust) usigma($cngr) vsigma() model(tre) marginal itera(100) rescale tolerance(1e-5) nrtolerance(1e-5) 

local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,jlms
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(NGR_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "NGR_tl9.dta", replace
restore
}
*
****************************************
*TZN
****************************************
{
preserve
keep if country=="TZN"
duplicates tag $idd, gen(tag)
drop if tag==0
tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
sfpanel lnout c.lnA##c.lnA c.lnL##c.lnL c.lnA#c.lnL $othinp tt2, distribution(exp) vce(robust) usigma($ctzn) vsigma() model(tre) marginal itera(100) rescale tolerance(1e-5) nrtolerance(1e-5) 

local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,jlms
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(TZN_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "TZN_tl9.dta", replace
restore
}
*
****************************************
*UGD
****************************************
{
preserve
keep if country=="UGD"
duplicates tag $idd, gen(tag)
drop if tag==0
tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
sfpanel lnout c.lnA##c.lnA c.lnL##c.lnL c.lnA#c.lnL $othinp tt2 tt3 tt4, distribution(exp) vce(robust) usigma($cugd) vsigma() model(tre) marginal itera(100) rescale tolerance(1e-5) nrtolerance(1e-5) 

local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,jlms
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(UGD_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "UGD_tl9.dta", replace
restore
}
*
****************************************
*GHANA
****************************************
{
preserve
keep if country=="GHN"

tab t, gen(tt)
xtset $idd t 
timer clear 1
timer on 1
frontier lnout ///
    c.lnA##c.lnA ///
    c.lnL##c.lnL ///
    c.lnA#c.lnL ///
    $othinp tt2 tt3, ///
    dist(exponential) ///
    vce(robust) ///
    iter(100) ///
    difficult uhet($cghn)
local conv=e(converged)
local ll=e(ll) 
local ite=e(iterations)
summ lnA
local mlan=r(mean)
summ lnL
local mlab=r(mean)
testnl (_b[lnA] + _b[lnL] ///
        + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
        + 1*_b[c.lnA#c.lnL]*(`mlan'+`mlab') = 1)	   
local chi2=r(chi2)
local pv=r(p)
predict inef,u
predict eff,te
summ eff
local TE=r(mean)
summ inef 
local TI=r(mean)
estat ic
mat k=r(S)
local akaike=k[1,5]
display `akaike'
timer off 1
timer list 1
local time=r(t1)
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append addtext(Chi2,`chi2',PValue,`pv',Technical_Efficiency,`TE',Technical_Inefficiency,`TI',Converged,`conv', Log simulated likelihood,`ll',Number of iterations,`ite', Akaike's information criterion,`akaike', seconds, `time') cttop(GHN_tl)
nlcom ///
 (eA: _b[lnA] + 2*_b[c.lnA#c.lnA]*`mlan' + _b[c.lnA#c.lnL]*`mlab')  ///
 (eL: _b[lnL] + 2*_b[c.lnL#c.lnL]*`mlab' + _b[c.lnA#c.lnL]*`mlan') ///
 (RTS: _b[lnA] + _b[lnL] ///
       + 2*_b[c.lnA#c.lnA]*`mlan' + 2*_b[c.lnL#c.lnL]*`mlab' ///
       + _b[c.lnA#c.lnL]*(`mlan'+`mlab')), post
outreg2 using "C:\Users\maur3\Dropbox\PC (2)\Desktop\Analisis_v9", excel bdec(3) append  
keep $idd t eff inef
save "GHN_tl9.dta", replace
restore
}

************************************************************************
* SUMMARY STATS - Table 4
************************************************************************
*Append datasets
{
use panel_1013Final.dta , clear
keep if country=="ETH"
merge 1:1 idd t using  "ETH_tl9.dta", keep(match)
tempfile t1
save `t1'

use panel_1013Final.dta , clear
keep if country=="MLW"
merge 1:1 idd t using  "MLW_tl9.dta", keep(match)
tempfile t2
save `t2'

use panel_1013Final.dta , clear
keep if country=="NGR"
merge 1:1 idd t using  "NGR_tl9.dta", keep(match)
tempfile t3
save `t3'

use panel_1013Final.dta , clear
keep if country=="TZN"
merge 1:1 idd t using  "TZN_tl9.dta", keep(match)
tempfile t4
save `t4'

use panel_1013Final.dta , clear
keep if country=="UGD"
merge 1:1 idd t using  "UGD_tl9.dta", keep(match)
tempfile t5
save `t5'

use panel_1013Final.dta , clear
keep if country=="GHN"
merge 1:1 idd t using  "GHN_tl9.dta", keep(match)
tempfile t6
save `t6'
}

clear
foreach num of numlist 1/6 {
	append using `t`num''	
}
bys country: summ output lnout farmsize lnA laborT lnL $othinp $ceth, sep(0)
tabstat output lnout farmsize lnA laborT lnL $othinp $ceth, by(country) c(s)
tabstat output lnout farmsize lnA laborT lnL $othinp $ceth, c(s)


************************************************************************
* TABLE AND PLOTS TECHNICAL EFFICIENCY - Table 7, Figures 4 and 5
************************************************************************
*Append datasets
{
use panel_1013R.dta , clear
keep if country=="ETH"
merge 1:1 idd t using  "ETH_tl9.dta", keep(match)
tempfile t1
save `t1'

use panel_1013R.dta , clear
keep if country=="MLW"
merge 1:1 idd t using  "MLW_tl9.dta", keep(match)
tempfile t2
save `t2'

use panel_1013R.dta , clear
keep if country=="NGR"
merge 1:1 idd t using  "NGR_tl9.dta", keep(match)
tempfile t3
save `t3'

use panel_1013R.dta , clear
keep if country=="TZN"
merge 1:1 idd t using  "TZN_tl9.dta", keep(match)
tempfile t4
save `t4'

use panel_1013R.dta , clear
keep if country=="UGD"
merge 1:1 idd t using  "UGD_tl9.dta", keep(match)
tempfile t5
save `t5'

use panel_1013R.dta , clear
keep if country=="GHN"
merge 1:1 idd t using  "GHN_tl9.dta", keep(match)
tempfile t6
save `t6'
}
clear
foreach num of numlist 1/6 {
	append using `t`num''	
}
{
replace eff=1 if eff>1 & eff!=.
gen country2=1 if country=="ETH"
replace country2=2 if country=="MLW"
replace country2=3 if country=="NGR"
replace country2=4 if country=="TZN"
replace country2=5 if country=="UGD"
replace country2=6 if country=="GHN"
label def country2 1 "Ethiopia" 2 "Malawi" 3 "Nigeria" 4 "Tanzania" 5 "Uganda" 6 "Ghana"
label val country2 country2
}
*Table TE summary stats (Table 7)
tabstat eff, by(country2) s(mean median sd min max skewness kurtosis) f(%6.2f)
*Plot TE by country (Fig4)
twoway (histogram eff, by(country2) bcolor(navy%80)xlabel(0(0.25)1)), scheme(white_brbg) 
graph export "TE_bycountry.png", replace 
*Plot TE by crop (Fig5)
preserve
bys country main_crop_std: egen obscrops=count(t)
keep if obscrops>200 
drop if main_crop_std==27
twoway (histogram eff, by(main_crop_std) bcolor(navy%80)xlabel(0(0.25)1)), scheme(white_brbg) 
restore
graph export "TE_bycrop.png", replace 



************************************************************************
* Quantile Effects - Figure 6
************************************************************************
*Append datasets
{
use panel_1013R.dta , clear
keep if country=="ETH"
merge 1:1 idd t using  "ETH_tl9.dta", keep(match)
tempfile t1
save `t1'

use panel_1013R.dta , clear
keep if country=="MLW"
merge 1:1 idd t using  "MLW_tl9.dta", keep(match)
tempfile t2
save `t2'

use panel_1013R.dta , clear
keep if country=="NGR"
merge 1:1 idd t using  "NGR_tl9.dta", keep(match)
tempfile t3
save `t3'

use panel_1013R.dta , clear
keep if country=="TZN"
merge 1:1 idd t using  "TZN_tl9.dta", keep(match)
tempfile t4
save `t4'

use panel_1013R.dta , clear
keep if country=="UGD"
merge 1:1 idd t using  "UGD_tl9.dta", keep(match)
tempfile t5
save `t5'

use panel_1013R.dta , clear
keep if country=="GHN"
merge 1:1 idd t using  "GHN_tl9.dta", keep(match)
tempfile t6
save `t6'
}
clear
foreach num of numlist 1/6 {
	append using `t`num''	
}
replace eff=. if eff>1 
gen ln_eff=ln(eff)

*Loop
global fir "prop10"
global sec "traderdens"
{
*------------------------------------------------------------
* 0) Create interactions prop10 x country and traderdens x country
*    (run once; can be done inside a do-file before the loop)
*------------------------------------------------------------
levelsof country, local(ctry)

foreach c of local ctry {
    gen first_`c'   = $fir * (country=="`c'")
    gen secon_`c'    = $sec * (country=="`c'")
    label var first_`c' "Commercialization x `c'"
    label var secon_`c'  "Trader density x `c'"
}
	
*----------------------------
* Settings
*----------------------------
local Q ".10 .25 .50 .75 .90"
tempname posth
tempfile qcoef

cap postclose `posth'
postfile `posth' str3 country q ///
    b_prop se_prop b_trad se_trad using `qcoef', replace

*----------------------------
* Loop over quantiles
*----------------------------
foreach tau of local Q {

     qreg   ln_eff ///
        first_* secon_* ///
        lnA lnL pesticides fertilizer ///
        agehead female ac_headeduc hhsize  ///
        i.year i.countryid, quantile(`tau') iterate(10000) 

    * Store country-specific slopes
    foreach c of local ctry {
        local b1  = _b[first_`c']
        local se1 = _se[first_`c']
        local b2  = _b[secon_`c']
        local se2 = _se[secon_`c']

        post `posth' ("`c'") (`tau') (`b1') (`se1') (`b2') (`se2')
    }
}

postclose `posth'
use `qcoef', clear

	
*----------------------------
* 95% CI bands (normal approx)
*----------------------------
gen ub_prop = b_prop + 1.96*se_prop
gen lb_prop = b_prop - 1.96*se_prop

gen ub_trad = b_trad + 1.96*se_trad
gen lb_trad = b_trad - 1.96*se_trad
}

*Plot
{
gen country2=1 if country=="ETH"
replace country2=2 if country=="MLW"
replace country2=3 if country=="NGR"
replace country2=4 if country=="TZN"
replace country2=5 if country=="UGD"
replace country2=6 if country=="GHN"

label def country2 1 "Ethiopia" 2 "Malawi" 3 "Nigeria" 4 "Tanzania" 5 "Uganda" 6 "Ghana"
label val country2 country2
}
twoway ///
    (rarea lb_prop ub_prop q, sort color(navy%15)) ///
    (line  b_prop q, sort lcolor(navy) lwidth(medthick)) ///
    (rarea lb_trad ub_trad q, sort color(maroon%15)) ///
    (line  b_trad q, sort lcolor(maroon) lwidth(medthick) lpattern(dash)) ///
    , ///
    xtitle("Quantile of log(TE)") ///
    ytitle("Coefficient") ///
    ylab(-0.25(.25)1.25,format(%5.2f) valuelabel) xlab(0(.25)1,format(%5.2f)) ///
    by(country2, cols(3) note("") title("Distributional effects on technical efficiency", margin(small))) ///
	yline(0, lcolor(gs6) lpattern(shortdash) lwidth(thin)) ///
    legend(order(2 "Commercialization" 4 "Local Trader density") ///
           position(12) ring(0) cols(2) size(*1.5)) ///
    scheme(tab2)
graph export "QuantilePlot_BYCOUNTRY_interactionsv9_gha.png", replace
